Bayesian Estimation for Land Surface Temperature 
Retrieval: The Nuisance of Emissivities 
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Abstract — An approach to the remote sensing of land 
surface temperature is developed using the methods of 
Bayesian inference. The starting point is the maximum en- 
tropy estimate for the posterior distribution of radiance in 
multiple bands. In order to convert this quantity to an es- 
timator for surface temperature and emissivity with Bayes' 
theorem, it is necessary to obtain the joint prior probability 
for surface temperature and emissivity, given available prior 
knowledge. The requirement that any pair of distinct ob- 
servers be able to relate their descriptions of radiance under 
arbitrary Lorentz transformations uniquely determines the 
prior probability. Perhaps surprisingly, surface temperature 
acts as a scale parameter, while emissivity acts as a location 
parameter, giving the prior probability 
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Given this result, it is a simple matter to construct estima- 
tors for surface temperature and emissivity. A Monte Carlo 
simulation of land surface temperature retrieval in selected 
MODIS bands is presented as an example of the utility of 
the approach. 

Keywords — Remote Sensing, Land Surface Temperature, 
Sea Surface Temperature. 



I. Introduction 

THIS paper derives the joint prior probability for sur- 
face temperature and emissivity for the land sur- 
face temperature retrieval problem in remote sensing. It 
presents analysis necessary for formulating a Bayesian ap- 
proach to that problem, together with a Monte Carlo simu- 
lation of land surface temperature (LST) and surface emis- 
sivity extractions. 

After a brief description of the problem and the method 
of attack, the maximum entropy estimator for the mis- 
match between sensed and forward model radiance is given. 
Next, the joint prior probability for surface temperature 
and emissivity is obtained. This quantity is required in 
order to construct a useable estimator for surface temper- 
ature and emissivity. With the prior probability in hand, 
it is a simple matter to construct expressions for the ex- 
pected values of surface temperature and emissivity con- 
sistent with sensor aperture radiances and available prior 
knowledge. Finally, a sample tcmperature-emissivity sepa- 
ration is presented using MODTRAN calculations both for 
the forward model and for simulated sensor radiances. 
©2004 The Aerospace Corporation 



II. The Temperature-Emissivity Separation 
Problem amd its Discontents 

Increasingly capable remote sensing technology has 
sparked interest in exploiting thermal emission from sur- 
faces both for remote sensing of surface temperature and 
of emissivity. On the one hand, surface temperature studies 
form a significant portion of the science goals of MODIS, 
ASTER, and MTI, while AVHRR has been used opera- 
tionally for sea surface and land surface temperature stud- 
ies for many years. On the other, the use of emissivity 
information in thermal portions of the spectrum for geo- 
logical remote sensing has grown rapidly in recent years, 
and is as prominent in the science goals of MODIS and 
similar instruments as is surface temperature. 

Accordingly, the problem of temperature-emissivity sep- 
aration merits close examination. The entirety of the in- 
formation about a radiating surface comes from the ther- 
mal radiation it emits, conventionally parameterized as the 
product of blackbody radiation at the surface temperature 
T and the emissivity e^, 



Is — ekBk{T) 
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at each wavenumber k. Suppose one observes a radiat- 
ing patch of a surface at each of n wavenumber intervals, 
and that one knows how to correct for the effects of line- 
of-sight (LOS) attenuation and contamination by radiance 
from other sources. One then has n equations of the form 
(0, but n+\ unknowns, including the surface temperature. 
In the absence of knowledge about T or from extraneous 
sources, one has an underdetermined problem. 

A variety of methods has been proposed for handling the 
temperature-emissivity separation (TES) problem|lj. In 
most approaches to this problem, simultaneous LST and 
band emissivity retrieval depends upon specifying an emis- 
sivity value in one or more reference bands. The MODIS 
Land Surface Temperature (LST) algorithmic seeks a 
pair of reference channels in a part of the thermal spec- 
trum in which the emissivity of natural surfaces displays 
very limited variation, and may therefore be regarded as 
known with good confidence. Multiband emissvities in- 
ferred on this basis are called "relative" emissivities |3_. Al- 
gorithms using this approach include the reference channel 
method^] and emissivity normalizationjS]. In the former, 
a value of emissivity is assumed for one band, and in the 
latter, an approximate surface temperature is obtained by 
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noting that cmissivity cannot exceed unity, in order to close 
the system of equations. Other relative emissivity retrieval 
methods include the temperature-independent spectral in- 
dex method|ni,[Z| and spectral ratios|H|- The study by Li 
et al. shows that all of these relative emissivity retrieval 
algorithms are closely related, and argues that they may 
be expected to show comparable performance. A different 
approach has been proposed for analysis of Multispectral 
Thermal Imager (MTI) dataj^I , in which radiances are col- 
lected from a surface with looks at nadir and 60 degrees 
off-nadir, assuming a known angular dependence of emis- 
sivity, in order to balance equations and unknowns. The 
generalized split-window LST algorithm^D] likewise uses 
dual looks in a regression-law based approach. The ba- 
sis of the "grey body emissivity" approach [TT] is the slow 
variation of emissivity with wavelength for certain natu- 
ral targets. The physics-based MODIS LST algorithmic 
exploits observations taken at day and at night, on the as- 
sumption that band emissivites do not change over periods 
of a few weeks. 

It is clear that the methods described depend upon a 
priori assumptions about the variation of emissivity, ei- 
ther with wavelength, or with look angle, or over time, 
from which one would like to be free. The work described 
in this paper uses Bayesian inference to retrieve estimates 
of surface parameters. This approach allows one to treat 
emissivities as "nuisance" parameters which may be inte- 
grated out of a posterior distribution function between par- 
simoniously chosen, and hence "uninformative," limits. 

It might appear odd to use as an approach to the sepa- 
ration of temperature and emissivity a Bayesian estimator 
which, in essence, allows one to ignore the actual value of 
emissivity. Equation shows that thermal radiance is 
linear in emissivity. However, the Planck function goes as 
a fairly high power of the temperature in the LWIR, and 
is close to exponential in temperature in the MWIR. Any 
roughness in the treatment of sensed radiance-as in allow- 
ing the assumed emissivity in the estimator to take on a 
wide range of values-may, therefore, be expected to lead to 
comparatively small errors in the inferred surface temper- 
ature. In fact, it turns out that the posterior distribution 
for surface temperature to be developed gives sharp limits 
to allowable surface temperature even in the presence of 
considerable latitude in the value of possible emissivities. 
In most cases, only a narrow range of surface temperature 
is consistent with the sensed radiance in multiple bands, 
whatever be assumed about emissivity. 

Once a reasonably good estimate of surface temperature 
is in hand, it is a simple matter to insert it into estimators 
for the individual band emissivities, and for the uncertainty 
in those values consistent with available knowledge. The 
a priori limits on emissivity may then be contracted, and 
a new estimate of surface temperature obtained. The ex- 
pectation values of surface temperature and emissivity may 
thus be refined iteratively. It is also possible to search for a 
(local) maximum for the posterior likelihood for these pa- 
rameters. Because the TES problem is underdetermined. 



this will not give a unique global maximum, but, given 
the insensitivity of surface temperature estimates to small 
emissivity errors, the local maximum may be expected to 
give results close to the physical values for the parameters 
of interest. 

III. Maximum Entropy Estimators for Surface 
Parameters 

Consider the problem of estimating surface temperature 
and emissivity from radiance detected by a remote sen- 
sor. The sensor supplies measurements of radiance / at 
the aperture. A forward model is required to compute the 
value of aperture radiance as a function of, among other 
things, the parameters we wish to extract. Assume ini- 
tially (for simplicity) that the sensor has fine spectral res- 
olution. The forward model radiance may be described at 
each wavenumber fc by a form of the Duntley equation 

Ipik) - ekBk{T)exp{^^) + P^Fj:{Q)exp{-^±) + iI{t, ^i) 

(2) 

ll (r, fi) and (0) are the upwelling diffuse radiance at 
nadir optical depth r (top of the atmosphere, or TOA, for 
spaceborne sensors; /i is the cosine of the angle with respect 
to zenith) and the downwelling irradiance at the surface, 
respectively. -Bfc(T) is the Planck function at surface tem- 
perature T. The emissivity is Cfc, and the surface reflectance 
Pk = 1— Efc. The form of (0) is what one would get assuming 
a Lambertian surface obeying Khirchoff's law. The anal- 
ysis presented below does not depend upon a Lambertian 
approximation to surface reflectances; in fact, it makes no 
assumption regarding their angular behavior |14| . In what 
follows it will be assumed that the only unknown quanti- 
ties in the preceding equation are T and e^. Generalization 
of the analysis which follows to the case of reflectance not 
equal to one minus the emissivity poses no difficulties. 

An estimator for the probability that, given observed 
radiances /(fc), the surface parameters T and take 
on particular values, is constructed in the following way 
dliCniidli {vide, also a related discussion in Landau 
and LifschitzfTH). pp. 343-5). The posterior probability for 
T and is given by Bayes' theorem as 

PiT,e, I /,if) = PiT,e, \ K)WljIhIl (3) 

where K denotes available knowledge. The quantity P{I \ 
K), the prior probability of the radiance /, may be ab- 
sorbed into an overall normalization and does not concern 
us further. It may be that the surface T is of interest, 
whatever the value of emissivity. In this case, one is free 
to denigrate as a "nuisance" parameter and integrate it 
out of 0, as will be done below. 

Consider the remaining factors in in turn, starting 
with the direct probability P{I \ T,ek, K) of observing ra- 
diance / given T, e^, and other a priori knowledge K. 
With aid of the forward model, it is possible to recast this 
quantity in more tractable form. By hypothesis, 

I{k) = lF{k) + eu (4) 
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where the error in spectral radiance is attributed to 
noise processes. The prior probabihty for the noise P{ek \ 
T, efe, K) is now obtained by a maximum entropy argument. 
If the noise power is assumed known, the noise probabihty 
is the function which maximizes the information-theoretic 
entropy subject to constraints imposed by the value of the 
noise power and overall normalization of probability, 
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P(e I K)log{P{e \ K))de - 

) 

/-t-oo r+oo 
e^P{e I K)de - A2 / P(e | K)de. 

The function maximizing ((^l is a Gaussian, 

P{e I K) - 
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(5) 



(6) 



where the Lagrange multipliers for noise power and nor- 
malization have been written in terms of the RMS noise 
radiance tr. Upon substituting for the noise term, the 
probability of detecting a radiance / given T, e, and noise 
a becomes 



P{I I T, e, a) — exp 
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dl 



(7) 



This is also the likelihood function for T and e. In order to 
formulate an estimator for T and e, it remains to find the 
prior probability 



P{T,e\K)^ f{T,e)dTde 



(8) 



(The appropriate prior for noise is known to be the Jeffreys 
form, but is omitted here because it is assumed that the 
noise contribution is known.) 

We now follow Jaynes' prescription for finding an unin- 
formative prior probability^ , [20] , |2] . Assume two equiv- 
alent observers record the same sensor aperture radiance 
originating as thermal radiation from a surface, and inter- 
pret it in terms of Planckian emission characterized by a 
surface temperature and emissivities, subject to LOS at- 
tenuation. Vladimir detects surface thermal emission I in 
a solid angle fi, and describes the surface with parameters 



T and e, and the attenuation with optical depth ^: 
I ^ ekBkiT)exp{~-) 



(9) 



He assigns prior probability in light of his knowledge re- 
garding the problem 



f{T,e)dTde 



(10) 



On the other hand, Estragon agrees with Vladimir on the 
definition of the Planck function, emissivity, and LOS at- 
tenuation, but describes the same situation with surface 
emission /' in il', and parameters T', e', and jp-, reporting 



I' = ei,B,,iT')exp{^^) 



(11) 



and assigning the prior probability 
g{T', e')dT'de' 



(12) 



In order for the pair to agree as to the form of the estimator, 
the priors must be related by 



g{T',e')dT'de' = J^^f{T, ejdTde 



where 



J — det 



d{T',e') 



(13) 



(14) 



a(r,e) 

is the Jacobian determinant for the transformation between 
descriptions in the parameter space. 

Assuming both are sober, Vladimir and Estragon must 
always be able to relate their descriptions of the sensed ra- 
diance by a Lorentz transformation. Let us consider active 
transformations for concreteness. Suppose that Vladimir 
wishes to describe events in a frame of reference moving 
at velocity (3 = v/c along the observation axis, denoted 
x, with respect to the frame preferred by Estragon. (It 
is convenient, although not actually necessary, to suppose 
that the x axis is also the axis of photon propagation.) 
Lorentz invariance requires that Vladimir's (unprimed) and 
Estragon's (primed) description of events be invariant un- 
der the Lorentz transformation given by 



(15) 
(16) 

, (17) 

The four-momentum of a photon travelling along the x-axis 



where 



x' — j{x — (3ct) 
t' ^-f{t- Px/c) 

1 
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/ hk/c \ 
hk 



V 



(18) 



Applying (|15|l and (|16|l to the components of (|18|l . we see 
that Estragon and Vladimir relate their description of fre- 
quency or wavenumber by 



fc' =7(1 -/3)fc 
= rjk. 



(19) 



How does the pair relate their descriptions of radiance? 
Let a bundle of Sn photons with mean energy p'^ — hck 
and uncertainty Sp^ — hcSk originate in a small area 6 A of 
the radiating surface in a small time interval St, collimated 
within a small solid angle 6fl, and propagate unattenu- 
ated to an observer. The surviving photons arriving at the 
observer's location comprise a collisionless photon gas. A 
single photon in the bundle occupies a phase space volume 



T4 X Vfc = SA{cSt) X h^k^SkSn 



(20) 



while the bundle occupies a 6 (5n-dimensional phase space 
volume 

Vphase = [V^Vk]'" . (21) 
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Equation (|21|l is invariant at any point on a photon trajec- 
tory. According to a standard result in statistical physics, 
Liouville's theorem (JH, pp. 9-10; 178), it has the same 
value at every point on that trajectory. So long as the 
photons remain coUisionless they can neither leave their 
original volume of phase space, nor enter another. How- 
ever Vladimir or Estragon choose to describe the patch of 
emitting surface SA, the time interval St, the solid angle 
interval SQ or the photon wavenumber k, they must agree 
as to the number of photons Sn in the bundle. Hence, both 
l|2()|l and the ratio 



N 



5n 



Sn 



v^Vk n^sAstk^Sksn 



(22) 



are invariant along a photon trajectory. 
Spectral radiance is defined as 



h 



d{energy) 



d{time)d{area)d{frequency)d{solidangle) 
Rewriting the radiance as 



hkSr 



SAStSkSn 



gives 



h . 

— ^ = invariant 
k-^ 



(23) 



(24) 



(25) 



for any component of the total radiance along a given line 
of sight 1221 ■ Equation H25|) has the same value in any frame 
of reference |23|. with two consequences for this problem: 
1. The Planck function obeys 



jy3 



BkiT) 



as may also be seen by direct substitution in 
1 fc3 



BkiT) 



exp 



hck 



ksT 



1 



(26) 



(27) 



2. Vladimir and Estragon must agree that the attenuated 
surface emission obeys 



fc'3 fc3 

Consider first the case of no attenuation, r = 0. Then 
e',,Bk'{r) _ eMT) 



(28) 



fc'3 



k^ 



or 



t't,,Bk.{T')^rfekBk{T). 
One also has, from H2t)|l . 

B^k{T') = rfBk{T'/f^). 

Combining and ||2U gives 

e',,Bk{T'/f^) = eMT). 



(29) 
(30) 

(31) 
(32) 



Now, Vladimir and Estragon also agree that, while they 
must lie between and 1, emissivities are otherwise com- 
pletely arbitrary functions of wavenumber, and by hypoth- 
esis have no dependence upon temperature. In 



BkiT) 



(33) 



"rjk 



the right-hand side can have no dependence upon T or T', 
while the left-hand side cannot be an arbitrary function of 
1] OT k. The only remaining possibility is 



Bk{T'/Ti) _ 
BkiT) 



const. 



(34) 



^rik 



The set of Lorentz transformations forms a group so 
this relation holds for the identity with (3 = 0, j = t] = 1. 
The constant must therefore equal unity. 
Next allow r to differ from zero in H28() . In 



k^e'^,Bk'iT') _ expi~ 
k'^ekBkiT) " 



-) 



expi~-^) 



(35) 



the left-hand side has, by hypothesis, no dependence upon 
LOS transmission, while the right-hand has no dependence 
upon surface properties so, again, both sides equal a con- 
stant, and, from (|26|l and 134|) . we find 



k%kBk'iT') ^ e',,B,kivT) 
k'^CkBkiT) ekV^BkiT) 



1 



or 



expi' 



expi- 



(36) 



(37) 



LOS attenuation does not affect the validity of H34II . 

Thus, the most general relation which respects a Lorentz 
transformation carrying wavenumber k to k' = rjk is 



T' = T]T 



The Jacobian is therefore 



(38) 
(39) 

(40) 



and 



fiT,ek)=vgi7^T,e'^,). (41) 

Invocation of the principle of indifference^ to assert Es- 
tragon and Vladimir must use the identical description of 
events, and thus assign the same prior probabilities. 



/(T,e)=g(r,e) 
leads to the functional equation 

fiT,ek) = vfi7^T,ek). 



(42) 



(43) 

^As given by Javnes IT^ . p. 128, in an extension of the original 
concept introduced by Laplace to encompass indifference between de- 
scriptions by distinct but equally cogent observers. 



5 



The solution of ijlSIl is 



const. 
T 



yielding 



fiT,e)dTde=^-^dTde 



(44) 



(45) 



for the prior probability. 

One now argues this form of the prior is least informative 
as to emissivity. No functional dependence upon a param- 
eter should enter the form of the prior probability that is 
not imposed by the requirements of invariance and indiffer- 
ence. Any such dependence would amount to the admission 
that we possess additional knowledge about emissivity be- 
yond that assumed. That is, H45(l is the unique choice of 
prior probability that assumes nothing about the value of 
ek beyond what is dictated by the problem statement. A 
standard argument (found, for example, in 15 , pp. 9-15) 
then shows that prior knowledge about limits on the value 
of emissivity should appear in the limits of integration used 
in constructing marginal distributions for T. 

Surface temperature thus obeys the Jeffreys prior, while 
emissivity obeys the Bayes prior. Both results may appear 
somewhat surprising, especially that for emissivity. From 
the manner in which it appears in the expression for radi- 
ance, one's naive expectation might be that emissivity is 
a scale parameter. However, the relation between the de- 
scription of emissivity as seen by Vladimir and Estragon 
more resembles what one would expect of a location pa- 
rameter: They must agree on the value of emissivity, but 
are free to assign it to different wavenumbers. 

The result just obtained will now be extended to the sit- 
uation in which radiance is sensed in bands wide enough 
that that it cannot be regarded as a function of wavenum- 
ber, but must be treated as an integral over a passband. 
One then writes, for the contribution of surface emission 
to the total radiance at the sensor aperture in band i. 



k2 rk2 

ekBk{T)exp{ — -)dk = ei / Bk{T)exp(- — )dk (46) 



It is always possible to do this by the mean value theo- 
rem for integrals, and it is frequently the case that the 
right-hand side of H4t)|l expresses all available knowledge 
concerning the radiant properties of the emitting surface. 
Vladimir describes the surface emission by 

Q / ' Bk{T)exp{-^)dk = f \kBk{T)exp{- — )dk (47) 
Jki M Jki M 

with 

0, fc < fci 
e-k = { £i,ki < k < k2 
0, fc > fc2 

while, by (|37|) - (|39|l . Estragon describes things by 



7^2 



e',,Bk'hT)expi-^)dk' 



7^2 



7fci 



Bk.{^T)exp{~^)dk' (48) 



with 



0, fc' < 7fci 
4' = { ^iilki < k' < 7fc2 
0, fc' > 7fc2 

Comparison of the two expressions for surface emission, 
(I47|l and (|48|l . leads to the immediate conclusion that the 
Jacobian connecting the two descriptions of surface tem- 
perature and band emissivity is 



J = det 



d{T\e') 
diT,e) 



1, 



and 



/(T,e)dm = ^rfT& 



(49) 



(50) 



once more. 

The result just obtained allows us to derive estimators 
for surface temperature and emissivity. The starting point 
is a calculation of the marginal posterior probability for T 
given observed radiance in a finite number of bands when 
the surface emissivity in band i is known to lie between 
£min{'i) and emax{i)- This quantity is computed for each 
band by integrating (j^J between these limits, upon insert- 
ing Q and (|45|) . Evaluating the integral requires complet- 
ing the square in the exponent of ([TJ. To accomplish this, 
define auxilliary quantities a, b, and c, obtained from jSJ: 



k2 
ki 



Bk{T)^-Fl{G)]exp{^ 
n 



)dk 



b = 6162 



with 

bi 



fci 



fei 



and 



Bk{T)--Ft{0)] exp{-^)dk 
) 



-Fl{^)exp{-^) + / I (T,Ai) \dk-h 

TT /i ' 



ii^,^(0)exp(-^)+4(r,M))dfc 



(51) 
(52) 

(53) 

(54) 
(55) 



Then (dropping subscript i for the moment) (O and Q 
give 

(ae2 



P(r,e|/,a)cx 



1 



27ro- 



exp 



be - 



2cr2 



(56) 



The marginal distribution obtained by integrating over the 
nuisance parameter e is 



P{T I /,cr) cx 



1 



27rcr 



exp 



be + c) 



2(t2 



de (57) 
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Completing the square in the exponent allows this to be 
written as 



P(T I /, cr) oc —=exp 

Ja 



[c - 6V4a] 



2a2 



^ i^^max : ^min) (58) 




for each band i. In (|59|) the error function is 

2 r 

erf{x) — —= \ exp (—t^) dt 
Vt^ Jq 



(59) 



(60) 



The joint posterior probability for observing radiances 
— I, n is 



^(r|/M=i,n,a) = []p(r|/„a). 



(61) 



Assuming T is known to lie between a minimum and a 
maximum, an estimator for T given radiance in band i is 



(r) = 



g-;rp(r|/„a)f^ 
/J;:;p(T|/„a)f 



(62) 



while a joint estimator for T given radiances in all n bands 
is 



dT 



(T) 



1 dT 



An estimator for the emissivity in band i is given by 

re:::n{T).e\h,a)d, 



(63) 



(64) 



This form has the advantage that estimates of the surface 
T are significantly less sensitive to discrepancies between 
sensed and modeled radiances than are estimates of emis- 
sivity. An estimate of T obtained from lf35|l - (|^ with un- 
informative limits on emissivity may be close enough in 
practice for accurate emissivity retrievals by H64|l . (Equa- 
tion (|MI) may be evaluated in closed form with elementary 
functions; however, the resulting expression is quite cum- 
bersome and is omitted here.) 

IV. Monte Carlo Simulation of MODIS Land 
Surface Temperature Retrieval 

A land surface temperature retrieval algorithm has been 
developed using the results just derived. While the intent 
of the work reported in this paper is to unshackle LST 
estimation from emissivity knowledge, the algorithm also 
retrieves emissivity estimates, and may be thought of as a 
TES algorithm if desired. It is intended to illustrate the 



application of Bayesian analysis to thermal remote sensing, 
and is purported to be optimal neither in execution speed 
nor accuracy. The requirement for only one forward model 
calculation per retrieval suggests that it will not impose an 
extreme computational burden in practice (even though the 
forward model to be described requires two MODTRAN 
calculations). The algorithm is used to simulate LST re- 
trieval from a notional cxoatmospheric sensor that records 
radiance from a patch of the Earth's surface at a specified 
signal-to-noise ratio(SNR). It is assumed that the domi- 
nant noise contribution arises from the shot noise of the 
radiance incident upon the sensor aperture. 

In outline, the algorithm works as follows. Equation 
gives the distribution of surface temperature consis- 
tent with observed radiances and the initial range of emis- 
sivities. This distribution typically differs from zero only 
within a narrow range about the true surface temperature. 
Within this range, H63|) is used for each of n bands and 
for all bands jointly to compute n -I- 1 separate estimates 
of the surface temperature. The actual surface tempera- 
ture is assumed to lie between the extreme values of this 
set of expected values, which now determine the allowable 
range. The joint temperature distribution and the various 
expectation values for surface temperature are next refined 
using the contracted range of a-priori credible surface tem- 
peratures in (|58|) . now calculated with a finer temperature 
mesh. After a few iterations of this procedure, the differ- 
ent surface temperature expectation values obtained from 
l|62|l and 163|l reliably converge to a single value lying close 
to the true surface temperature. A convergence radius rj 
for the different estimates of 0.01 K was used. Emissivities 
are then obtained by substitution of the joint surface tem- 
perature estimate into the expression for band emissivity 
expectation values, H64|l . Equation H56|l being a Gaussian 
distribution, it is possible to refine the a-priori limits on 
credible emissivites by specifying a threshold of m stan- 
dard deviations. Six is used for the examples presented 
here. The revised a-priori emissivity limits and surface 
temperature limits may then be used to restart the entire 
sequence just outlined, if desired. This additional itera- 
tive loop was repeated once in the simulation presented 
here. Surface temperature and emissivity values show only 
marginal changes as a result of the second iteration, indi- 
cating convergence of the retrieval. 

The starting point is the posterior distribution for sur- 
face temperature To compute it, the coefficients a(T), 
6(r), and c are obtained as a function of surface temper- 
ature. Isaacs two-stream M0DTRAN4 calculations with 
SALB=1.0 supply those forward- model spectral quantities 
independent of surface temperature: Attenuation along the 
line-of-sight, downwelling radiance at the surface, and up- 
welling radiance at TOA. The surface thermal emission 
component is computed directly from the Planck function 
and spectral attenuation. A further MODTRAN calcula- 
tion with SALB=0.0 is used to obtain estimates of scat- 
tered solar radiance and of the "scattered thermal" radi- 
ance compoment of the total radiance returned by MOD- 



7 



TRAN. 

The calculation of a(T) and b{T) used in the retrievals 
departs from (|51|I - H53|I in one regard. Approximation of 
the total TOA radiance computed by MODTRAN with 
the (never exact) Duntley equation becomes increasingly 
inaccurate as the surface temperature increases. The dom- 
inant contribution to the discrepancy is the surface emis- 
sion portion of the scattered thermal radiance. Subtract- 
ing an estimate of this term gives a corrected TOA ra- 
diance in much better accord with the predictions of the 
Duntley equation. The estimate is calculated as a function 
of the unknown surface emissivity and (necessarily erro- 
neous) forward model boundary temperature, using the dif- 
ference between the scattered thermal contributions com- 
puted with SALB=1.0 and SALB=0.0. Rather than sub- 
tracting the approximate scattered surface thermal radi- 
ance from the total MODTRAN radiance, the correction 
was implemented in a mathematically equivalent way by 
adding that portion of the scattered thermal radiance com- 
ponent linear in surface emissivity to the surface thermal 
emission terms in H51|l - H53|l . 

It should be noted that while the forward model as- 
sumes knowledge of atmospheric parameters, MODTRAN- 
computed quantities used in the forward model have no de- 
pendence upon true surface temperature or emissivity. The 
boundary temperature parameter used in MODTRAN af- 
fects the forward radiance calculations primarily through 
the scattered thermal contribution. In addition, MOD- 
TRAN adjusts the atmospheric temperature profile in the 
lower zones to interpolate smoothly between surface con- 
ditions and a fiducial layer in the atmosphere. For these 
reasons the forward model boundary temperature should 
not differ greatly from a physically reasonable value. 

The algorithm is executed according to these steps: 

1. Perform the forward model radiative transfer calcula- 
tions with MODTRAN. 

2. Calculate the individual band posterior probabilities, 
and the joint posterior probability over all bands, as a func- 
tion of surface temperature with H62|l and H63(l , contracting 
as necessary the range of T to that giving nonvanishing 
joint posterior probability in H61|) . 

3. Calculate expectation values for surface temperature 
over the posterior probabiltiy for each band individually, 
and over the joint posterior probabiltity. This calculation 
gives n + 1 surface temperature estimates. 

4. Perform convergence test max\{Ti) — {Tj)\ < rj over all 
pairs of surface temperature estimates. If the convergence 
test is satisfied, proceed to step 5. Otherwise, iterate by 
repeating steps 2-4. 

5. Compute expectation values for band emissivities us- 
ing 

6. Adjust emim £max to ±m Standard deviations about 
(fi) for each band. 

(7. Repeat steps 2-6 if desired.) 

Monte Carlo simulations of LST retrieval in a selected 
subset of MODIS bands illustrate the performance of the 
algorithm. The bands chosen appear in Table 1. MOD- 



TRAN calculations are used both as simulated TOA radi- 
ances in MODIS bands and as the forward model. Each 
Monte Carlo realization of TOA radiance is calculated us- 
ing a mid-latitude summer atmosphere with MODTRAN 
parameters listed in Table 2 selected as uniform deviates 
within the limits shown, including "true" surface T and 
band emissivities. It is unlikely that the atmospheric pro- 
file or other parameters required to specify a formard radia- 
tive transfer model will be reliably known to high accuracy. 
In order to simulate the effect of imperfect knowledge on 
the forward model, a second draw of random numbers is 
used to introduce errors in the fallible forward model as 
shown in the last column. Thus, for examples, the surface 
visibility, which MODTRAN uses to parameterize aerosol 
effects, is chosen to lie between 5 and 30 km for each Monte 
Carlo realization. A random error of up to ±4 km is added 
to this value of visibility for use in calculation of the fal- 
lible forward model for each realization. The MODTRAN 
model default water vapor profile, given as (grams precip- 
itable water)/ (kilograms air), is randomly scaled between 
the limits shown (subject to the constraint that relative 
humidity cannot exceed 100%). The range of perturbed 
forward model parameters is truncated at the limits spec- 
ified in Table 2. The forward model cannot, of course, in- 
corporate knowledge of the true surface temperature, but 
it does require an initial guess for that quantity. This guess 
is obtained by varying the forward model boundary tem- 
peraature randomly from "truth" by ±20^, without tnm- 
cation. Inclusion of the scattered thermal radiance con- 
tribution in the forward model notably improves retrieval 
accuracy, despite boundary temperature uncertaintiies of 
this magnitude in the forward model. 

Table 1. MODIS bands used in simulations 

MODIS band wavelength limits notional snr 

20 3.660-3.840 fi 350 

22 3.929-3.989 350 

23 4.020-4.080 350 
29 8.400-8.700 /i 1000 

31 10.870-11.280 A* 1000 

32 11.770-12.270 ii 1000 



Table 2. Monte Carlo parameter ranges 



MODTRAN 


minimum 


maximum 


fallible forward 


parameter 


value 


value 


model error 


nadir view angle 


125 dcg 


180 deg 


±0.125 deg 


surface visibility 


5 km 


30 km 


± 4 km 


column water vapor 


0.33 X MLS 


1.00 


±0.2 


thin cirrus altitude 


8 km 


12 km 


±0.5 km 


thin cirrus thickness 


1 m 


20 m 


±25 m 


thin cirrus opacity 


0.05/fcm 


0.2/A;m 


±0.025/fem 


solar azimuth 


deg 


90 deg 


±0.125 deg 


solar elevation 


20 deg 


60 deg 


±0.125 deg 


viewng azimuth 


deg 


90 deg 


±0.125 deg 


viewing elevation 


35 deg 


90 deg 


±0.125 deg 


surface T 


268 deg K 


328 deg K 


N/A 
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Next, the simulated MODIS TOA radiances are contam- 
inated with notional sensor noise simulated as a zero-mean 
Gaussian random process with standard deviation equal to 
the noise equivalent radiance NEAR. The algorithm also 
requires an estimated variance for the noise radiance (a 
in (jnj), which should be of order NEAR. NEAR is pa- 
rameterized in terms of a signal-to-noise ratio. SNR was 
chosen to lie on the low end of values inferred from MODIS 
NEAT values pS] with the aid of the following estimate. 
The error in TOA radiance from noise sources is estimated 
as, roughly, 



Table 3. Monte Carlo simulation results: mean errors, 
standard deviations 



5I^—5T + 0{5Tf 



(65) 



dT 



NEAT 



Tk . 



(66) 



where 

Is = el Bk{T)exp{~ — )dk 

JAk M 

is the attenuated surface thermal emission at TOA, leading 
to 



SNR- 



61 



f dlogjls) 
[ dT 



NEAT] 



(67) 



This quantity was computed for each Monte Carlo realiza- 
tion; the SNR values used for the retrievals, which apppear 
in Table 1 , conservatively underestimate 1)67(1 . Performance 
of the algorithm appears not too sensitive to the exact noise 
contamination added to the simulated band radiances, nor 
to the exact noise variance assumed in the retrieval, as long 
as neither is grossly erroneous. In fact, it is possible to ad- 
just the assumed value of a in the estimator to contract 
or expand the range of viable surface temperatures consis- 
tent with sensor radiances without disastrously biasing the 
retrieval, as described below. 

One thousand Monte Carlo realizations each were calcu- 
lated for day and night, with a mid-latitude summer at- 
mosphere. Generous bounds for the initial a-priori limits 
on LST and band emissivity were assumed, subject to the 
physical upper bound on surface emissivity: 



200K <T < 500K 
0.75 < e, < 0.99 



(68) 



Note that both limits for LST lie outside the range sampled 
by the Monte Carlo draws. 

Mean errors and error standard deviations for the re- 
trieved surface temperatures and band emissivities, with 
respect to "true" values, appear in Table 3. In the majority 
of cases, acceptable estimates of LST and band emissvities 
were obtained using all six bands from Table 2, with the 
SNR chosen to equal the assumed noise variance in the sim- 
ulations which appears in Table 1. However, in about 4% of 
the simulations it proved impossible to find an acceptable 
solution with all six bands in this manner. The solution 
instead tended to badly erroneous values (e.g., T < lOOK, 
and emissivities pegged at the limits of the prior). Inspec- 
tion of the posterior probabilities from H62|l revealed that in 



Case 



Day 



Night 



LST error(K) 
£20 error 
€22 error 
£23 error 
£29 error 
€31 error 
£32 error 



-0.25 ± 1.23 
-0.004 ± 0.022 
-0.009 ± 0.034 
-0.008 ± 0.048 
-0.004 ±0.031 
-0.005 ± 0.023 
-0.007 ± 0.028 



-0.31 ± 1.11 
-0.003 ± 0.035 
-1-0.001 ± 0.034 
-0.007 ± 0.038 
-0.003 ± 0.022 
-0.005 ± 0.022 
-0.006 ±0.029 



these anomalous cases one or more of the individual band 
posterior probabilities fails to overlap significantly with the 
product of the remaining posterior distributions, leading 
to a joint posterior probability that effectively vanishes. A 
number of remedies is available when this difficulty arises. 
The number of successful retrievals rises sharply when the 
range of the band emissivity prior is expanded to 0.7-0.999, 
but at the cost of reducing their overall accuracy somewhat. 
Experimentation shows that, in all of the anomalous cases, 
it is possible to get a satisfactory LST retrieval with some 
three-band subset of the original six. The LST so obtained 
can then be inserted into the expectation value for band 
emissivity to yield emisssivities for all six bands. Finally, 
the support of the joint posterior probability can be broad- 
ened by increasing the noise radiance a in its calculation. 
This last approach was used to obtain Table 3. The effect 
on retrievals of increasing a for subsets of the bands dif- 
fered negligibly from that of increasing a for all bands by 
the same factor. The most intractable of the anomalous 
cases (one each nighttime and daytime) required increas- 
ing cr by a factor of 7.0 in order to obtain an acceptable 
solution. 

Examination of the statistic for retrieval errors shows 
that the estimator is (slightly) biased for normal, and sig- 
nificantly biased for anomalous, retrievals. For the 963 
normal nighttime retrievals the mean and stardard devi- 
ation of the surface temperature error are ST = — 0.26 ± 
1.06K, with = 1-06 per degree of freedom. For 958 
normal daytime retrievals, the corresponding figures are 
ST = -0.21 ± 1.18K and =1-03 per degree of free- 
dom. The values for 37 anomalous nighttime retrievals are 
ST — —1.71 ± 1.32K and — 2.64 per degree of freedom, 
with ST = -1.08± 1.92K and =1-29 per degree of free- 
dom, for 42 anomalous daytime cases. However, the distri- 
bution of errors is very accurately Gaussian (as should be 
expected). If the mean error is subtracted before calculat- 
ing x^, the result is identically 0.999 per degree of freedom 
for surface temperature (and all six band emissivities) for 
all retrievals. 

Forward models for the anomalous cases systematically 
have large errors in one or more of the randomly-varied 
simulation parameters. Thus, the surface visibility and the 
boundary temperature are both more likely to lie near the 
limits of their range than in the middle. In particular, the 
column water vapor scaling factor is over three times as 
likely to exceed 1.1 for anomalous retrievals as for normal 
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ones (25/37 vs 203/963 nighttime; 27/42 vs 182/958 day- 
time), with only three daytime and two nighttime anoma- 
lous retrievals occurring for a column water vapor scaling 
less than unity. The frequency of anomalous retrievals ap- 
pears, to some extent, to be an artifact of inserting large 
errors in the forward model to approximately simulate im- 
perfect knowledge of atmospheric conditions. 

In any event, all 2000 Monte Carlo realizations led to a 
successsful retrieval of both LST and six band emissivities. 

V. Discussion 

Points which should be addressed in further develop- 
ments of practical algorithms: 

1. It appears that this approach to TES works largely 
because the range of plausible surface temperature values 
consistent with band radiances and an uninformative range 
of band emissivity is quite constricted, as a consequence of 
the strong temperature dependence of the Planck function. 
It turns out not to be terribly difficult to get a tempera- 
ture estimate that is close enough to truth that it can be 
inserted into the least sophisticated imaginable estimator 
for band emissivity (|64|l . and still lead to acceptably accu- 
rate results. 

Once the algorithm has gotten to an iteration in which 
the current range of temperature and band emissivities is 
restricted to a neighborhood sufficiently close to the true 
values that the posterior distribution is jointly Gaussian 
in T as well as in the e^, it is apparent both that con- 
vergence to the true values will occur as assumed in this 
algorithm, and that these values will maximize the likeli- 
hood. However, at present, there is no proof in hand that 
the procedure outlined above actually converges, or that, 
given that it does, it converges to the true surface temper- 
ature and emissivity combination. It appears to do both 
to good accuracy in practice. 

However, the algorithm did-initially-fail to converge to 
an acceptable solution in about 4% of the realizations. As 
recounted in the previous section, it proved possible in ev- 
ery case to adapt the search strategy so as to successfully 
retrieve both LST and emissivities for all bands. The suc- 
cessful recovery strategies all had the effect of maximizing 
the numerical joint posterior probability, by some combi- 
nation of 1) eliminating from the estimator band poste- 
rior probabilities whose effective nonvanishing support did 
not intersect that of the joint probability of the remaining 
bands, or forcing intersection by broadening the support 
of the outlier posterior probabilities by 2) loosening limits 
on the prior, or 3) increasing the noise radiance parameter 
assumed in the retrieval. 

The solutions obtained are, in any event, not unique, 
because the TES problem is underdetermined. Considered 
as a surrogate for maximum likelihood solutions, the al- 
gorithm solutions approximate only local maxima, and it 
might be possible to find maxima which give very poor 
account of temperature and emissivity. This has not hap- 
pened in simulations performed to date. 

2. The algorithm as presently formulated appears to be 



unnecessarily complicated. It seems certain that its opera- 
tion can be significantly streamlined. For practical applica- 
tions, it will be necessary to eliminate redundant elements 
of the calculation. 

3. The model for band radiances in this memo treated 
them independently, apart from the prior knowledge that 
the surface temperature for all bands must be the same. 
Bretthorst ^SliEliEl 1^^^ addressed problems involving 
more sophisticated models for observations, in an approach 
which would appear to offer real advantages in the present 
context. 

4. Perhaps the least satisfactory feature of this algo- 
rithmic approach is its dependence upon an accurate for- 
ward radiance model. To the extent MODTRAN can be 
regarded as supplying radiance estimates which are zero- 
mean error estimates of the true radiance, the effect of 
radiance prediction error on this algorithm may simply be 
regarded as a contribution to the noise variance. But in real 
life, a forward model can be expected to have systematic 
errors that need not originate as unbiased stationary Gaus- 
sian processes. The question whch has been addressed in 
this work is: Given an accurate forward model (in the sense 
just described), what surface temperature and band emis- 
sivities are consistent with observed radiances and knowl- 
edge of their error statistics? A harder question, which will 
be the focus of further developments, is: Given a fallible, 
but reasonably accurate, forward model, what surface tem- 
peratures and band emissivities can possibly be consistent 
with observations and available knowledge, no matter what 
the forward model error, so long as it falls within known 
limits? 

VI. Conclusion 

A simple argument, based on inherent physical symme- 
tries that the description of surface thermal emission must 
obey, leads to the appropriate prior probabilities for surface 
temperature and emissivity. These lead to the maximum 
entropy estimator for the mismatch between sensed and 
modeled radiance in the presence of noise, from which an 
estimator of surface temperature may be constructed that 
treats emissivity as a nuisance parameter. MODTRAN- 
based simulations show that temperature-emissivity sepa- 
ration is successfully performed by iteration between the 
temperature estimator, and a similar estimator for surface 
emissivity. 
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